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Abstract 

A  hybrid  computational  method  for  solving  boundary-value  problems  is 
introduced  which  combines  features  of  the  meshless  hp-cloud  methods  with 
features  of  conventional  finite  elements.  The  method  admits  straightforward 
nonuniform  hp-type  approximations,  easy  implementation  of  essential  bound¬ 
ary  conditions,  is  robust  under  severe  distortions  of  the  mesh,  and  can  deliver 
exponential  rates  of  convergence.  Results  of  numerical  experiments  are  pre¬ 
sented. 


1  The  Method 

Recent  developments  made  in  the  context  of  meshless  methods  have  demonstrated 
the  simplicity  of  adding  hierarchical  refinements  to  a  low  order  set  of  shape  functions 
which  satisfy  the  partition  of  unity  (PU)  requirement  [1—4, 6,  7].  In  particular,  in 
the  hp-cloud  method  introduced  in  [2],  one  covers  the  domain  Q,  of  the  solution  of  a 
boundary- value  problem  with  a  collection  of  open  sets  C  U”_iu;c),  the  sets 

uja  being  the  clouds,  and  constructs  on  the  clouds  a  set  of  global  basis  functions  (fa 
which  form  a  PU  on  fl: 

a=l 
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One  can  then  build  spectral  (p-type)  approxinaations  by  constructing  products  of 
the  Kpa  with  higher-order  spectral  approximations  to  produce  enriched  basis  functions 
(poLp,  Lp  being  a  polynomial  of  degree  p.  This  idea  has  been  used  successfully  in  gen¬ 
erating  exponentially  convergent  approximations  of  elliptic  boundary-value  problems 
in  which  convergence  is  obtained  by  appropriate  /i-refinement,  h  =  maxoha,  ha  being 
the  diameters  of  the  clouds  cUq,  or  p-enrichment,  p  being  the  order  of  the  polynomial 
carried  by  the  basis  functions  (paLp  [2-4, 6, 7].  Here,  the  basic  building  block  for  such 
approximation  is  the  PU,  and,  as  observed  in  [1],  this  partition  of  unity  can  be  fur¬ 
nished  by  a  conventional  finite  element  method.  While  such  a  PU  on  a  finite  element 
mesh  will  destroy  the  “meshless”  quality  of  the  approach,  ample  compensation  for 
this  loses  is  provided  by  a  number  of  advantages  over  conventional  /ip-finite  element 
methods. 

Consider,  for  example,  the  conventional  finite  element  meshes  of  triangles  or 
quadrilaterals  shown  in  Figs.  1  and  2,  respectively,  on  which  continuous  global  ba¬ 
sis  functions  (shape  functions)  Na  are  constructed  at  each  nodal  point  Xa,  o.  = 
1, 2, . . . ,  n.  These  functions  are  such  that 

n 

^  Na{x)  =  1,  at  any  x  E  Cl 

a=\ 

and  thus  form  a  PU.  By  setting  pa  =  we  can  build  /ip-cloud  approximations 
on  this  PU  and  thereby  produce  a  cloud  basis  that  has  all  of  the  useful  convergence 
properties  of  the  /ip-clouds  but  which  combines  features  of  conventional  FEM’s,  such 
as  exhibiting  the  Kronecker-delta  property  at  boundary  nodes.  Some  examples  of 
cloud-type  basis  functions  are  shown  in  Figs.  3,  4  and  5.  Notice  from  the  figures  that 
(see  also  Figs.  6,  7,  8) 

•  The  clouds  uia  need  not  be  disks  or  even  convex  polygons, 

•  The  mesh  parameter  h  (or  ha)  is  the  diameter  of  the  cloud  (the  patch  surround¬ 
ing  node  Xa)  and  not  the  diameter  of  an  element  in  the  cloud 

Figure  8  shows  a  conventional  hierarchic  field  as  used  today  in  many  codes 
and  introduced  in  the  manner  described  originally  by  Zienkiewicz,  et  al  [11]  and  later 
elaborated  on  by  others  [8,10]. 

(a)  The  conventional  hierarchic  polynomial  are  introduced  using  local,  element 
based,  co-ordinates.  These  provide  complete  Cartesian  polynomials  only  up  to  the 
linear  terms  when  isoparametric  distortion  is  introduced.  This  completeness  can  be 
extended  to  quadrilateral  elements  using  quadratic  terms  only  by  introduction  of  the 
9  node  expansion  (as  shown  in  [12]).  For  higher  order  terms,  performance  of  elements 
may  well  deteriorate;  viz  [5]. 
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Figure  5:  Higher  order  hierarchical  shape  function  built  from  the  product  of  the 
bilinear  shape  function  shown  in  Fig.  3  and  the  monomial 


Figure  6:  Overlapping  patchs  corresponding  to  clouds  uja  and  ujp.  Polynomials  of 
differing  degree  pa  and  pp  can  be  associated  with  nodes  at  Xa  and  xp  so  as  to  produce 
non-uniform  Ap-cloud/finite  element  approximations. 


Linear 


Linear  +  Quadratic 
Linear  +  Quadratic  +  Cubic 


Figure  7:  Cloud-based  hierarchic  fields  in  two  dimensional  space 


Quadratic 
Quadratic  +  Cubic 

Figure  8:  Conventional  hierarchic  fields  in  two  dimensional  space 


On  the  contrary,  the  new  cloud-bcised  hierarchic  forms  give  the  terms  of 
quadratic  and  higher  order  expansion  in  terms  of  Cartesian  coordinates  through¬ 
out  and  will  always  retain  the  accuracy  corresponding  to  the  spectral  order.  Indeed 
some  saving  in  degrees  of  freedom  necessary  for  a  given  accuracy  is  available. 

(b)  The  new  hierarchic  form  concentrates  all  the  unknown  degrees  of  freedom 
at  corner  nodes  of  the  elements  (cf.  Fig.  7).  This  ensures  a  more  compact  band 
structure  than  that  arising  from  the  conventional  hierarchic  form. 

(c)  Remarkably,  the  structure  of  this  approximation  allows  the  use  of  different 
values  of  the  spectral  order  p  on  each  cloud.  Thus,  other  that  basic  connectivity  of  the 
low-order  FEM  mesh,  none  of  the  complications  of  high-order  constraint  conditions 
used  in  hp-  finite  element  methods  are  necessary.  Non-uniform  h  and  p  approximation 
can  be  easily  accommodated  over  the  mesh,  as  suggested  in  Fig.  6. 

(d)  Further,  the  elimination  of  corner  degrees  of  freedom  follows  precisely  the 
same  pattern  as  that  used  for  the  underlying  linear  element  repeating  the  same  elim¬ 
ination  steps  in  the  equation  solver.  This  means  that  vectorization/parallelization  of 
the  algorithms  is  relatively  easy.  What  is  more,  if  adaptive  steps  are  used  in  the  re¬ 
finement,  insertion  of  higher  order  terms  at  any  node  does  not  alter  the  computational 
sequence  and  can  be  simply  accomplished. 

(e)  The  new  process  presents  no  difficulties  on  the  boundaries  if  Dirichlet 
(prescribed  displacement)  boundaries  occur  and  no  higher  order  terms  are  introduced 
in  boundary  nodes  than  a  simple  specification  of  corner  displacement  suffices. 

With  higher  order  terms  on  boundary  nodes  and  with  locally  curved  bound¬ 
aries,  additional  constraints  may  be  required  (exactly  as  in  standard  hierarchic  ele¬ 
ments). 

For  Neumann  (prescribed  traction)  boundaries,  no  such  difficulties  arise  pro¬ 
viding  the  integration  is  carried  out  on  the  exact  curves. 

1.1  Construction  of  Linearly  Independent  Basis  Functions 

The  cloud-based  hp  finite  element  basis  functions  are  defined  above  as 

(pat  • —  NaZf,’ 

where  Na  is  a  finite  element  shape  function  and  Li  is  a  polynomial  of  degree  p.  Since 
the  shape  functions  form  a  PU 


=  E^-Li  =  LiY^Na  =  Li  (1) 

cx  a  a 

Therefore,  the  polynomials  Li  can  be  recovered  through  linear  combinations  of  the 
cloud  basis  functions  (pai-  The  finite  element  shape  functions  have  the  property  that 
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exist  Gxa,  Oyo,  €  R,  o  =  1, . . . ,  n  such  that 

^  y  ®xo  ~  ^ 

a 

y"!  Oyo  —  J/ 

a 

Now  if  we  take  L,  —  x  or  L,  =  y  we  have  from  (1)  that 

^  yVa^;  =  X  ^  NaV  =  y 

Or  Or 

Therefore  the  set  {A^„,  A^o3^}o=i  not  linearly  independent.  This  can  be  avoided 
by  a  careful  choice  of  the  functions  Li  used  to  build  the  cloud  basis  functions.  For 
example,  the  elements  from  the  space  span{Na}  must  not  be  used  to  build  the  cloud 
basis  functions. 


2  Numerical  Experiments 

2.1  Illustration  of  the  Basic  Ideas 

As  a  first  numerical  experiment,  we  consider  the  solution  of  the  following  boundary- 
value  problem 

Find  u  such  that 


—An  =  0  in  =  (— 1, 1)  x  (— 1, 1)  (2) 

where  x  +  iy  =  z  E  C  and  Re  denotes  the  real  part  of.  The  value  of  u  is  fixed  at 
(1,1)  in  order  to  make  the  solution  unique. 

Figure  9  depicts  the  solution  u.  The  initial  finite  element  discretization  is 
represented  in  Fig.  10.  It  is  composed  of  only  two  quadrilateral  bi-linear  finite 
elements.  The  corresponding  finite  element  solution  obtained  with  this  discretization 
is  depicted  in  Fig.  11  which  has  an  error  in  the  energy  norm  of 

=  79.06% 

iFlIf; 

Quadratic  shape  functions  like  those  depicted  in  Figs.  4  and  5  are  then  added 
to  node  0,  as  shown  in  Fig.  12.  The  corresponding  finite  element  solution  is  shown 
in  Fig.  13. 
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ure  10:  Initial  finite  element  mesh  composed  of  two  bi-linear  elements. 


Figure  13:  Finite  element  solution  obtained  using  the  discretization  of  Fig.  12. 


Figure  14  shows  the  finite  element  discretization  in  which  nodes  0, 1,3  and  4 
have  quadratic  approximation  associated  with  them,  while  nodes  2  and  5  have  bi¬ 
linear  approximations.  The  corresponding  finite  element  solution  is  shown  in  Fig.  15. 
The  discretization  error  associated  with  this  discretization  is 

7  =  34-78% 

II“IIe 


Figure  16  shows  the  finite  element  solution  obtained  using  quadratic  shape 
function  at  all  nodes  of  the  discretization.  As  expected,  the  error  for  this  discretization 
is  of 

11^  ~ 


w||j 


=  c?(io-*) 


It  is  emphasized  that  the  enrichment  of  the  finite  element  spaces,  as  described 
above,  is  done  on  a  nodal  basis  and  the  polynomial  order  associated  with  a  node  does 
not  depend  on  the  polynomial  order  associated  with  neighboring  nodes. 

The  insensibility  of  the  cloud-based  approximations  to  element  distortion  is 
next  demonstrated.  Figure  17  shows  the  finite  element  discretization  used.  All  nodes 
have  quadratic  cloud-based  approximations.  Figure  18  shows  the  corresponding  finite 
element  solution.  As  before,  we  get 

=  O(10-«) 
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Figure  16:  Finite  element  solution  obtained  using  quadratic  shape  functions  at  all 
nodes  of  the  mesh. 
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Figure  17:  Distorted  mesh  of  quadratic  cloud-finite  elements. 
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Figure  18:  Finite  element  solution  obtained  using  the  discretization  of  Fig.  17. 


2.2  hp  adaptivity 

W’e  shall  now  consider  the  problem  of  an  L-shaped  plane  elastic  body  loaded  by  the 
tractions  associated  with  the  following  stress  field 

t7„(r,0)  =  Ar^-*[(2-(5(AH-l))cos(A-l)^-(A-l)cos(A-3)0] 

=  Ar^-^[(2  +  (5(A  +  l))cos(A-l)0  +  (A-l)cos(A-3)0]  (4) 

a,y{r,6)  =  Ar^-‘[(A-l)sin(A-3)^  +  (5(A+l)sin(A-l)0] 

where  (r,  is  the  polar  coordinate  system  shown  in  Fig.  19,  A  =  0.544483  737, 
Q  =  0.543  075579. 

The  stress  field  (4)  corresponds  to  the  first  term  of  the  symmetric  part  of  the 
expansion  of  the  elasticity  solution  in  the  neighborhood  of  the  corner  A  shown  in  Fig. 
19  [9]. 

Plane  strain  conditions,  unity  thickness  and  Poisson’s  ratio  of  0.3  are  assumed. 
The  strain  energy  of  the  exact  solution  is  given  by  [9] 

a}^ 

S{u)  =  4.154  54423— 

E 

where  E  is  the  modulus  of  elasticity  and  a  is  the  dimension  shown  in  Fig.  19.  The 
values  E  =  a  =  1  are  eissumed  in  the  calculations. 
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Figure  19:  L-shaped  elastic  body. 

Two  sequences  of  discretizations,  Si  and  <^25  are  used  to  solve  this  problem.  In 
the  former,  the  uniform  mesh  shown  in  Fig.  20  is  used  and  the  polynomial  order  of 
the  approximations  ranges  from  1  to  8.  A  strongly  graded  mesh,  as  shown  in  Fig.  21, 
and  non-uniform  p  distributions  are  used  in  the  second  sequence  of  discretizations. 
Geometric  factors  q  =  0.10  and  q  =  0.15  for  the  size  of  the  elements  are  used  (cf.  Fig. 
21).  Figures  22  and  23  show  the  p  distribution  used  in  the  fourth  step  of  this  sequence. 
The  polynomial  order  of  the  clouds  decrease  linearly  towards  the  singularity  while  the 
size  of  the  finite  elements  decrease  geometrically. 

The  relative  error,  measured  in  the  energy  norm,  versus  the  number  of  degrees 
of  freedom  is  shown  in  Fig.  24  for  the  sequences  of  discretizations  Si  and  S2  (with 
q  =  0.10  and  q  =  0.15).  As  expected,  the  uniform  mesh  gives  an  algebraic  rate  of 
convergence  while  the  strongly  graded  meshes  lead  to  an  exponential  rate  of  conver¬ 
gence.  This  kind  of  behavior  is  typical  of  hp-  finite  element  methods.  However,  the 
construction  of  non-uniform  h-  and  p-  discretizations  in  a  cloud  based  framework  is 
considerably  more  straightforward  than  in  conventional  hp-  finite  element  methods. 
For  this  problem,  the  strongly  graded  mesh  with  q  =  0.10  gives  slightly  better  results 
than  the  case  g  =  0.15 
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Figure  22:  Polynomial  order  associated  with  the  clouds  at  the  fourth  step  of  the 
sequence  of  discretizations  82- 


Figure  23:  Zoom  at  clouds  near  the  re-entrant  corner  of  Fig.  22. 
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L-Shaped  Domain:  Cloud-FEM  Results 


Number  of  degrees  of  freedom 

Figure  24:  Convergence  in  the  energy  norm  for  uniform  and  non-uniform  cloud-based 
hp-  distributions. 
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3  Conclusions 


In  this  investigation,  we  have  explored  the  use  of  lower-order  finite  element  approx¬ 
imations  to  generate  partitions  of  unity  on  which  hierarchical  hp-  cloud  approxima¬ 
tions  can  be  constructed.  The  resulting  methodology  has  a  number  of  useful  features. 
Among  these  are  that  non-uniform  hp-  meshing  with  variable  and  hierarchical  order  p 
over  clouds  can  be  easily  generated.  The  spectral  convergence  of  p-  and  hp-  methods 
is  retained  and  the  method  is  very  robust,  the  accuracy  being  quite  insensitive  to 
mesh  distortion.  Also,  by  building  cloud  approximations  on  finite  element  meshes, 
Dirichlet  boundary  conditions  are  easily  handled. 

The  hp-  convergence  properties  seem  to  differ  from  traditional  p-  version  ele¬ 
ments,  but  exponential  convergence  is  attained.  Applications  to  problems  with  singu¬ 
larities  are  easily  handled  using  cloud  schemes.  In  all,  this  hybrid  finite-element /cloud 
methodology  appears  to  have  a  number  of  useful  and  attractive  features  that  could 
prove  to  be  important  in  broad  engineering  applications. 
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